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ABSTRACT 

High order accurate centered flux approximations used in the computation of numerical 
solutions to nonlinear partial differential equations produce large oscillations in regions of 
sharp transitions. In this paper, we present a new class of filtering methods denoted by ENO- 
LS (Essentially Nonoscillatory Least Squares) which constructs an upgraded filtered solution 
that is close to the physically correct weak solution of the original evolution equation. Our 
method relies on the evaluation of a least squares polynomial approximation to oscillatory 
data using a set of points which is determined via the ENO framework. 

Numerical results are given in one and two space dimensions for both scalar and systems 
of hyperbolic conservation laws. Computational running time, efficiency and robustness of 
the method are illustrated in various examples such as Riemann initial data for both Burgers’ 
and Euler’s equations of gas dynamics. In all standard cases the filtered solution appears to 
converge numerically to the correct solution of the original problem. Some interesting results 
based on nonstandard central difference schemes, which exactly preserve entropy, and have 
been recently shown generally not to be weakly convergent to a solution of the conservation 
law, are also obtained using our filters. 
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1 Introduction 


Numerical improvements in the computation of high order accurate numerical solutions to 
nonlinear hyperbolic conservation laws have been recently obtained. Hence, following Total 
Variation Diminishing (TVD) schemes, the Essentially Nonoscillatory (ENO) method has 
been introduced and proved to be very efficient in the computation of high accurate numerical 
solutions for several types of physical problems including Computational Fluid Dynamics 
(CFD) problems or front propagation using the Hamilton-Jacobi framework. However, these 
higl^liccurate methods use a lot of computational time. For that reason, filtering methods 
were developed, beginning in the late eighties. The first one, described by B. Engquist 
and B. Sjogreen in [1], uses simple TVD and conservation properties to correct nonphysical 
spurious oscillations from one time step to another. The correction step consists in pushing 
numerical data points up or down to an acceptable level while preserving conservation. In [5], 
we presented a new class of filtering methods of any order of accuracy. Our method relies 
on switching fluxes at locations in which spurious oscillations are detected. This method 
was observed to be very efficient and its cost relatively low since high order ENO fluxes are 
evaluated at a few points only- a central difference method is used most often. 

In this paper, we investigate some interesting computational properties of centered schemes 
after numerical oscillations have developed and propagated for some time. We define a new 
class of filtering methods that can be applied to highly oscillatory numerical solutions. This 
relies on the construction of an ENO stencil of points ([2, 7, 8]) which is fitted with high de- 
gree polynomials from a least squares procedure. Our numerical filter is capable of smoothing 
oscillations having large amplitude and high frequency, but without removing sharp singu- 
larities which are crucial components of these solutions. Furthermore, the filtered solution 
seems to retain the oscillatory solution properties of the unfiltered schemes in some special 
“entropy preserving” cases as defined in [3] by J. Goodman and P. D. Lax, and in [6] by 
J. Liu and D. Levermore. We investigate some numerical examples using several centered 
approximations in order to illustrate the former property. The main conclusion indicates 
that for standard central differences, our filtered solution always converges accurately to the 
strong limit, whereas the predicted oscillatory behavior is retained even after using our filter 
in the examples of [3, 6]. 

Our main test problems will be inviscid Burgers’ equation and the inviscid Euler equations 
of gas dynamics. We first consider Burgers’ equation: 

u, + (^\ = o, (1) 

with smooth initial condition U(x, 0) = Uq(x), Uq <E C 00 ^, 1), and periodic boundary con- 
ditions. We will discuss the main properties of the numerical solution obtained from some 
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schemes based on the approximate fluxes: 


F m = + uj), ( 2 ) 

F m = sWh + U P) ~ hW* + £7-.). W 

*<+*= + + 555(^3 + ^). W 

= im „ + £/^, +1 + uj), ( 5 ) 

■^j+J ~ jt/jl/j+i • (6) 


The fluxes (2,3,4) are just standard central differencing of second, fourth and sixth order 
of accuracy, respectively; while (5,6) are the interesting examples analyzed in [6] and in [3], 
respectively. The oscillatory solution obtained from any of these fluxes is then corrected by 
the ENO-LS method, preserving the formal order of accuracy. 

The Euler equations of gas dynamics: 


U f + (F(U)) X = 0, (7) 

U(*,0) = Uo(z), 

are to be solved for t > 0 and x in some interval f) with appropriate boundary conditions, 
where 


F(U) = uU + (0 ,p,vp) T 

and U = (p,q,p), p is density, q is momentum, v is velocity, and p is the pressure. In this 
work, we use conventional second, fourth and sixth order central differencing with ENO-LS 
post processing applied to Euler’s equations. See [4] for an analysis of the oscillatory Von 
Neumann-Richtmyer scheme approximating Euler’s equation, [9]. 

2 The ENO-LS Filter, Algorithms 

The ENO-LS method mimics the construction of ENO polynomials but without involving 
the evolution equations. In short, we follow the algorithm just below: 

Algorithm 2.1 • 1.) Compute N times the numerical solution of 

Ut + A{V) X = 0 

U(x, 0 ) = U 0 (x) 

i.e we let V® = Uo(xj), and compute Vj N = I(V°,NAt), for all j = l,...,n, where I 
is the solution semigroup operator that transforms the initial pointwise data V.° to V) jV 
after N iteration time steps. 
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• 2.) Filter the numerical solution computed in step L) by an iteration procedure similar 
to Jacobi or Gauss-Seidel elliptic solvers; first let: 


W° = V- N 

3 3 

for all j = 0, n, then make use of primitive variables: 

u T+i = E w i" Ax " 

i=Q 

and finally construct a sequence } m =o ... m so that 

u ™+i = E(U m ,U m+1 ), 

where M is defined from, the stopping criterion: 

\\W m+1 - W m || < e, 


for m = 0 , M — 1, where e is a small parameter of order Ax a and 




C'fti-C'Jl i 
2 J 2 


A Xj 


Finally, let Vj v = Wf 1 , for j — 1, 
• 3.) Go to step 1.) unless t = t max . 


( 8 ) 

( 9 ) 


( 10 ) 


We notice that relation (10) and the use of primitive variables (8) implies the conservation 
property of the sequence {Wj}j=i,..., n , he the resulting finite difference scheme is always in 
conservation form. Moreover, the number of correction steps M can be initially fixed so that 
the ratio is as large as desired. The operator it 1 is a non trivial linear combination of 
, for some j , as in point Jacobi or Gauss Seidel method. In the Jacobi procedure 

we have: 


and either 







if j varies from 0 to n, or 
in the reverse direction; 


for the Gauss-Seidel method. 
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An important property of the ENO-LS algorithm comes from the fact that the corrected 
solution Vj N satisfies a conservation equation. To see that, we first use the relation (10) 
which can be rewriten as: 

(E»h(U H , rt} ...,U iH+ , t )- U? +i ) - (B'-iWi-i-r. ~ t£i) 


w m +l _ iym 
J ) 


A Xj 


( 11 ) 

and then discuss the construction of the least squares process. In (11), the pair of integers 
(r-j-jS-t) limit the width of the stencil used in the evaluation of the least squares polynomial. 
The appropriate stencil is defined as in the ENO algorithms (refer to [7], [8]). We briefly 
indicate the main steps leading to such a polynomial: We first compute the divided differences 
table of W m and define the ENO stencil of points in the region which is the smoothest for 
successive space derivatives of W m . We denote the ENO stencil in the set (x ; _ r , ..., Xj+ 3 ), 
where r + s = p — 1, and p is the number of data points that we want to take into account 
in the evaluation of the least squares polynomial. Hence if we denote this polynomial by 

= ( 12 ) 

*= o 

where (<^ 0 , V^) are the basis functions of some polynomial space of degree q , then the 

unknown coefficients Y = are solutions of the linear system: 

C j+ * Y = F j+ *, 

where is a (q + 1) x (q + 1) square matrix, F J+ ^ is a q + 1 column vector, and both 

can be computed from the basis functions: 

Tli=-r l Pk( x j+l+i) i Pl( x j+±+i)‘> 

F k 2 = 72U-T Uj+^+iVkiXj+l+i)- 

The updated value U m V follows from letting x = x j+ i_ on the previously constructed LS 
polynomial: 

try = r + i ( x,. + p = *)• 


»'=0 


The global conservation feature follows provided that we write: 

uzv = t 


t—-r 


where / = 0 or 1 depending on whether the Jacobi or the Gauss-Seidel method is used, and 

a i+l+t = ^2^2F> itk ^ipi(Xj + i)(pk(xj + i +t ), 


«= o k=0 
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J+k 


where = (C' J+ *) _1 . 

Note first that the coefficients {c/TF } t= _ r * form a sequence of bounded real numbers, 

J "t" 2 * * 

and second, that the basis functions ...<p q ) can be appropriately chosen as a sequence of 
Chebyshev polynomials or some other set of orthogonal polynomials. The last choice permits 
us to compute directly the coefficients a J+ 2 without inverting the matrix C J+ 2 : 


j+ l ^Vi(x j+ L)<Pi(x j+ l +t ) 

a i+H = E 

1=0 


cif 1 


This yields a fast algorithm since matrix inversions are no longer needed. 

We now present an extension of algorithm 2.1 to two space dimensions. The simplest 
possible extension would be to apply the previous algorithm in two sweeps. The first one 
will freeze one coordinate and correct the oscillatory solution with respect to the other free 
variable. The second one will simply reverse the role of each variable. This method, while 
simple, has difficulties near curved shocks. We shall use instead a fully two dimensional ENO- 
LS filtering algorithm. The latter will provide the construction of least squares polynomials 
in two space dimensions using a set of points which is chosen as the intersection of one 
dimensional ENO stencils of data points in each separate direction. The algorithm below 
describes our procedure. 


Algorithm 2.2 • 1.) Compute N times the numerical solution of 


U t + A(U) X + B(U) y = 0 

U(x,y, 0) = U Q (x, y) 

using a very simple numerical method and then let Vf- = Uo(xj,yi), and Vfj = 
I(V°,NAt), for all j = l,...,ni, and i = 1,...,U2, where I is the solution semigroup 
operator that we have already encountered in previous algorithm. 


2.) Let Wfj = Vfj , and filter M times the primitive variable 


hj 


= £ W^AxtAy,, 


l,k = 0 


by a Jacobi or point Gauss-Seidel iteration procedure , i.e perform: 


(13) 



E(U m ,U m+l ), 


for positive integers m. Iterate as long as ||W m+1 — W m \\ < e, m = 0,..., M — 1, and 
finally let the filtered solution: = Wfj , for j = l,...,nj, i = l,...,n 2 - 

The construction of the two dimensional least squares polynomial is as follows: 
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2.1 ) At the location (x i+ i,y j+ i), compute the ENO stencils of points 
{(•£»■— r*> J/j)j (®i+l+i*> Vj')} > (Tnd {(*£«> Vj — Tyfi •••» (•*-«> Vj + l+Sy )}> for r X T &X 1 Til 
and r y + s y + 1 = n 2 , where n i and n 2 are the preset maximum number of points along 
the x and y axis. These sets of points form two segments crossing at ( Xi,yj ). Then, 
define the x and y ENO stencils of points along these y and x segments, respectively. 
The two dimensional ENO stencil is taken as the intersection of the union of the pre- 
defined x and y one dimensional ENO stencils (refer to figure (1))- The least squares 
polynomial is then simply defined on this two dimensional ENO stencil and is set to 
pi'+bi+b^x^) = Ybk=oYkTk{x,y)- Again, the unknowns coefficients Yk, k = 


are computed by solving the linear system C^ +2 ’ J+2 ^Y = F^ +2,J+2 ^. 
S.S) Let 


2.3) Recover the conservative solution 


W . TO . +1 = 

h] 


vrX M + U ?_ tv , - u~V H - 


AxAy 


(14) 


for m = 0, ...,M — 1. 

• 3.) Go to step 1.) unless t = t maT 


An example of two dimensional ENO stencil is given in figure 1. The main interesting 
feature of such construction is based on the localization of the least oscillatory part of the 
solution within the much larger rectangle [xj- rx ,Xj +Sx ] x [y,-r v , ih'+aj- 

In section 4, we investigate the two dimensional Burgers’ equation and study numerical 
propagation of a shock along the radial axis. With centered fluxes, some spurious numerical 
oscillations propagate in the direction of the flow; however the two dimensional ENO-LS 
filtering method is able to filter all these oscillations while still giving the correct location of 
the curved shock wave. 

Next, we consider hyperbolic systems of conservation laws exemplified by the inviscid 
Eulers’ equations of compressible gas dynamics: 

Ut + F{U) X = 0, 

U(x, 0) = U 0 (x) 

for which there exists a complete set of real eigenvectors and eigenvalues; i.e VF(U) = 
P _1 AP, where A is a diagonal matrix with entries Ai < A 2 ... < Aj, and P is a matrix whose 
columns define a complete set of eigenvectors of the system. Note that the eigenvalues can be 
multiple, which is the case for the Euler equations of gas dynamics in two space dimensions. 


6 



Using the eigenvector decomposition, a field by field approach is used, i.e the fluxes will be 
corrected in each fields (our filter sometimes failed to remove oscillations near strong shocks 
when applied to the conserved quantities). Briefly, we proceed as follows: 


Algorithm 2.3 (Field by field ENO-LS method) • 2-1.) LetW m = (Wp, ..., W™) T , 
and define the primitive variable W™Ax. 

• 2-2.) Compute the left and right eigenvectors \ for k — 1 ,...,d (d is the 

number of equations) of the Roe matrix (see e.g. [2]) AfWf ,W^ +1 ) , and decompose 
the primitive vector V™ k along each characteristic field: 


m,(k) _ 



for which we have frozen the index j in order to get the same decomposition for all 
neighboring points Xj + i + /(j.p for l(k) — — r(k ), ..., +s(&) involved in the calculation of 
the least squares polynomial which is constructed in step 2-3.). 


• 2-3.) Select an ENO stencil of p points, i.e {xj_ r (fc), ...,Xj+ s (fc)}, for r(k) + s(k) + l = p; 
and then define the least squares polynomial P^ + 2 ) of degree q so that: 


g.m+l,(k) _ p(j+L), m,(k) 

j+l \ i+t- 




™>(k) ) 

i+j+»( fc )’ J+a'" 
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• 2-4-) Transform back the filtered field by field solution to the primitive vector: 


u 

-*7771 + 1 + 

v iU ~ a j+i K j+i > 


i+t t-s Hi 

2 k= 1 2 


and finally recover the desired physical variables: 

■ym+l _ -y-m+1 
wm+l _ £+ 2 1 

j Ax 

• 2-5,) Iterate until the stopping criterion 


( 15 ) 


||W m+1 — W m || < e, 

for m = 0 , M — 1, is reached; and finally let = W^. 

Note that this algorithm can be extended to two space dimensions by correcting separately 
oscillatory fields involved in the x and y fluxes, respectively. Moreover, this algorithm does 
not make use of the evolution equations but does use the eigenfunction decomposition in 
order to track efficiently the propagation of spurious oscillations. 

To conclude this section, we indicate that we mainly supposed that the numerical oscil- 
lations always propagate with the flow speed along local characteristic fields and that the 
amplitudes of such oscillations are not too large so that the oscillatory solution does not be- 
come unbounded, e.g., negative density and/or pressure is not allowed. In all our numerical 
experiments, we had to turn the filter on not only to recover an acceptable final solution but 
also to reduce the amplitude and frequency of spurious oscillations during the calculations. 
Hence, we usually preset the value of the ratio in the numerical experiments just after 
singularities have developed. 

3 Numerical Convergence Study 

In this section, we investigate the numerical convergence of the approximate solution com- 
puted via the ENO-LS algorithm given an initial oscillatory solution which has been eval- 
uated from one of the standard centered fluxes (2), (3), (4); or from one of the ’’entropy 
conserving” fluxes (5) or (6). As test problem, we study first the numerical evolution of 
Burgers’ equation (1) in one space dimension: 

U 2 

^ + (%-), = 0, 

with initial condition /7(ar,0) = sin 27ra:, in the domain [0,1], and extend the solution by 
periodicity outside 0 and 1. A shock wave develops at t — % at the point x = 
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We first display in figure (2) the solution using 100 cells and a CFL condition of | 
for a second order centered difference (CD) scheme using the flux (2), and then correct 
that solution by applying 7 iterations of the ENO-LS algorithm. Note first that the exact 
solution is perfectly recovered with exact location of the singularity because of the use 
of primitive variables in algorithm 2.1; second, these results are obtained for the set of 
coefficients (p, q) = (13,2), which implies that the corrected solution is only second order 
accurate in smooth transition areas; third, the upper right plot displays the corrected solution 
when both (3,2)CD and (3,2)ENO schemes are sequentially used. In that case, we correct by 
5 ENO iterations after every 25 CD steps. This method for filtering oscillations is also very 
efficient for one dimensional problems (refer to table 1); however, many more ENO iterations 
are needed in two dimensions. Basically, 10 ENO iterations are needed after every 10 CD 
steps in order to recover an accurate solution for the two dimensional Burgers’ example 4.3, 
which is a quite expensive technique compared to the overall cost induced by the ENO-LS 
method. 

We shall discuss the results obtained with larger values of q in section 4 in which we 
investigate numerical order of accuracy of the filtered solutions from the ENO-LS algorithm. 
Also, a similar study is investigated as the number of evaluation points p increases. 

In figure (3), we plot the numerical solution before and after the filter steps when 
n = 1000. Again, convergence to the physical solution is reached within 10 Gauss Sei- 
del iterations. Note that this number increases by a factor of nearly 2 when Jacobi iterations 
are used. 

In figures (4) and (5), we show the filtered solutions which were initially computed using 
the standard (3,4)CD and (3,6)CD schemes (fluxes (3) and 4)), respectively. Note that the 
corrected solution is well reconstructed in smooth regions while a smearing of the shock over 
about 10 cells is observed. This smearing can be primarily explained from the high (fourth) 
degree least squares polynomials ( q = 4) used in those experiments. 

Our second numerical test is devoted to showing that the filtered solution computed via 
the numerical flux (5): 

Fj+i = iWU + U > U M + UJ). 

does not converge to the physical solution. It was observed in [6] that this scheme preserves 
mass and entropy (which in this case is taken to be \U 2 ). It was also shown there that the 
numerical solution does not converge to a weak solution of the original problem. Basically, 
as the mesh size tends to zero, some spurious oscillations are still visible near the shock and 
cannot be removed. We ran the last two previous experiments but used the approximate flux 
(5). We plot the numerical results in the figures (6) and (7) when n = 100 and n = 1000, 
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respectively. For 100 cells (figure (6)), the oscillations near the shock are all smoothed out; 
however, the location of the shock is smeared over about 10 cells. This results were obtained 
for the same set of coefficients (p, q) = (13,2). Indeed, the solution is not well reconstructed. 
As the number of cells increases, the filtered solution is well reconstructed except near the 
shock. Numerical oscillations are still visible and cannot be removed. Note that the results 
visualized in figure (7) are given after 16 Gauss Seidel iterations which is much more than 
we needed when standard centered differences are taken. 

Our final Burgers’ equation test deals with a similar convergence failure property to the 
correct physical solution when the flux approximation (6) is implemented. The numerical 
flux is: 

This scheme again conserves both mass and entropy (this time the entropy is taken to be 
log Uj), and it was shown in [3] that the numerical solution does not converge to a weak 
limit of the original problem as the stepsize Ax tends to zero. In the numerical experiments, 
we noticed that the amplitude of the oscillations grew very fast and was rapidly becoming 
unbounded. The results are plotted in the figures (8) and (9) for n = 100 and n — 1000, 
respectively. Note that, in the case n = 100, the filtered solution is again not very well 
reconstructed, and for 1000 cells, some oscillations are still visible after 20 ENO-LS iterations 
on the right of the shock. Moreover these oscillations could not be removed, even after many 
additional filtering iterations. 

Thus our techniques have been observed to construct converging sequences of filtered 
numerical solutions towards the expected solution given initial oscillatory data that has 
been computed from a central difference flux approximation. Our hope is now to show that 
our method is not only robust but is also fast. This is the topic of last section. 

4 Time Efficiency of ENO-LS Algorithms 

In this section, we want to test the ENO-LS filtering method for several test problems 
involving nonlinear hyperbolic systems of conservation laws. We will focus our attention 
on comparing precisely the CPU times of the ENO-LS method versus more classical and 
filtered methods. Among them, we will consider the straightforward central difference (CD) 
method, the expensive ENO (Essentially Nonoscillatory) technique [7, 8], and our FM scheme 
(Filtering Method) [5], We will run three examples for ID and 2D Burgers’ equation, and 
for Eulers’ equations of compressible gas dynamics. 


12 











Type of Scheme 

Order of Accuracy 

CPU time xlO"* 

# of corrections, comments 

CD 


0.56 

CD = Centered Differences 

CD 

(3,4) 

0.72 

- 

ENO-RF 

(3,2) 

1.11 

RF = Roe Fix, Entropy fix 

ENO-RF 


1.68 

- 

FM 

(3,2) 

0.76 

FM = Filtering Method, 3% Corrections 

FM 

(3,4) 

1.07 

4% Corrections 

FM 

(3,2) 

1.40 

100% Corrections = full ENO 

FM 

(3,4) 

2.04 

100% Corrections 

CD+ENO (40,5,20) 


0.74 

40CD, 5ENO, 20 final ENO 


(3,2) 

0.58 

40CD, 1ENO, 10 final ENO 

CD+ENO (20,7,10) 

(3,4) 

0.99 

20CD, 7ENO, 10 final ENO 


(2) 

1.84 

( p } q) = (7, 2), LU Inversion 

ENO-LS (7,2) 

(2) 

0.65 

Orthogonal polynomials 

ENO-LS (7,3) 

(3) 

2.5 

LU Inversion 

ENO-LS (7,3) 

(3) 

0.89 

Orthogonal polynomials 


Table 1: CPU times of CD, ENO, FM and ENO-LS methods. 


4.1 ID Burgers’ 

We first compare the time efficiency of several numerical schemes for the ID Burgers’ equation 
(1). In table (1), we show the average computational time per iteration for 100 mesh points 
for the CD, ENO, and FM schemes with several correction angles (see [5]), and for various 
combinations, performing alternatively some centered difference and some ENO steps. The 
notation CD+ENO (40,5,20) simply means that the calculation starts with 40 CD steps, 
followed by 5 full ENO iterations, and back to the centered scheme; the last 20 steps of 
the calculations are finally performed by the ENO method in order to recover an acceptable 
solution. Note that all coefficients displayed in this table are tuned so that the numerical 
solution is high order accurate in smooth regions and no spurious oscillations are detected. 

The contents of this table need a few comments. First of all, the fastest algorithm is the 
one based on central differences. This is indeed not surprising since only one Fortran instruc- 
tion is needed in the coding of the approximate flux. Second, postprocessing a numerical 
solution computed from the CD scheme by an ENO method can be very efficient for lower 
orders. For a second order method, only one ENO iteration after every 40 CD steps has to be 
implemented in order to reduce sufficiently the amplitude of oscillations. Note however that 
for the fourth order method, 7 ENO iterations were needed every after only 20 CD steps. 
In fact, if these spurious oscillations are not regularly cut off, the final ENO iterations may 
not recover the correct numerical solution. Third, the ENO-LS filtering method is the most 
costly when full LU inversion of the C 3+ ? matrix is performed. However, when orthogonal 
basis functions are introduced, the ENO-LS method is competitive with respect to the fast 
CD scheme. Note moreover that ENO-LS correction steps have to be performed at a few 
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Type of Scheme 

(q.p) 

# of iterations 

L 1 and L°° orders. 

CD (3,2) 

(2,5) 

12 

1.81 and 1.63 

CD (3,2) 

(2,9) 

9 

1.75 and 1.50 

CD (3,2) 

(2,13) 

7 

2.19 and 1.9 

CD (3,4) 

(3,5) 

15 

1.59 and 1.40 

CD (3,4) 

(3,9) 

12 

2.41 and 2.33 

CD (3,4) 

(3,13) 

7 

2.80 and 2.79 


Table 2: Local L 1 and L°° order of accuracy. 


times only. Finally, after running many experiments, we noticed that if the ratio £ becomes 
too large, then shocks have a tendency to spread over a large number of cells. Again, there 
is a compromise that needs to be reached for fast convergence; this depends on the large 
value of the ratio and the approximation of the shape near shocks for which p needs to 
be slightly smaller. Note that in most of experiments, the ratio ^ £ [4,6] was optimal. 

Now, we want to measure the order of accuracy of the ENO-LS solution. To do so, we 
measure from computations the L 1 and L°° errors in the slabs [0.10,0.24] and [0.26,0.30]. 
Numerical orders are shown in table 2. 

Several comments about these results are now discussed. First of all, as the number of 
evaluation points increases, the better the quality of the results and the faster convergence 
is reached. Indeed, we have to pay the price of higher computations which are required to 
construct the coefficients involved in the matrix and in the vector F J+ ^ . On the other 

hand, faster calculations can be obtained provided that the values of p and q differ only 
slightly, i .e p = q + l, l = 1,2,3, ..., but local accuracy becomes obviously poor. Again, the 
optimal ratio seems to belong to the interval [4,6]. 

4.2 Example 2. Euler Equations of Gas Dynamics 

The second test problem is devoted to the Euler equations of gas dynamics in one space 
dimension. We consider the initial condition given in example 8 of [8], that is: 

p = 3.857143,9 = 2.629369, p = 10.3333333 when x < — 4 
p — 1 + e sin 5x, q = 0,p = 1. when x > — 4 

When £ — 0, a pure Mach 3 shock is moving to the right from the initial discontinuity 
x = 4. When e = 0.2, we not only have a Mach 3 shock propagating to the right, but have as 
well a succession of weaker rarefaction and shock waves propagating to the left. Numerical 
results for the ENO and FM methods of high order of accuracy can be found in [8] and in [5], 
respectively. We ran the same problem using 240 CD iterations and then plotted the results 
in figure (10). Next, we correct this highly oscillatory numerical solution by performing 7 
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Figure 10: (3,2)CD scheme. 

ENO-LS correction steps. In this experiment, we use (p, q) = (13,3). The filtered solution 
is visualized in figure (11). Note that the pressure, velocity, and entropy are quite well 
reconstructed, whereas the density is not perfectly recovered near the strong shock for which 
some physical oscillations should remain ( refer to [10, 11]). However, we should note the 
remarkable improvement obtained from the solutions displayed in figures (10) and (11). 

In figures (12) and (13), we use in sequence 50 CD steps and only one ENO iteration 
for second and fourth order methods before correcting the final oscillatory solution. The 
oscillatory solution is now postprocessed by 3 ENO-LS iterations. The filtered numerical 
results now look fairly similar to those shown in [8, 5]. 

4.3 Example 3. 2D Burgers’ Equation 

The last example is devoted to the two dimensional Burgers’ equation to be solved in the 
square domain [— 1, 1] x [—1, 1] with initial condition Uo{x,y) = sin27rr, where r = T y 2 , 
for r < and U 0 (x,y) = 0 outside the disc r = J. In figure (14), we visualize the solutions 
obtained by the (3,2)CD scheme, followed by 4 iterations of ENO-LS correction steps. In 
that experiment, we use (p xt p y ,q) = (6,6,2), so that the local rectangle in which the two 
dimensional ENO-LS stencil of points is taken contains 36 mesh points. Note that a speed up 
factor of nearly 1.5 is obtained by using the two dimensional ENO-LS method instead of the 
one dimensional splitting version. In this numerical example, the amplitude of the numerical 
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Figure 11: (3,2)CD scheme + 7 ENO-LS method (p,q)=(13,3). 

oscillations which propagate radially away from the initial shock location was approximately 
^ of the strength of the shock, yet the ENO-LS method recovered the nonoscillatory accurate 
solution quite well. 

5 Concluding Remarks 

The notion of using an essentially nonoscillatory (ENO) least squares filter to remove spu- 
rious oscillations of data generated by numerical overshoot appears promising. Tuning of 
parameters is, however, still required. Future work will deal with this issue. 
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-PRESSURE- 



- (3,2)CD+ENO (50,1) + 5 final (ENO-LS) (p,q)=(13,3) - 

Figure 12: (3,2)(CD,ENO) schemes + 3 ENO-LS method. 



- (3,4)CD+ENO (50,1) + 3 final (ENO-LS) (p,q)=(13,4) - 


Figure 13: (3,4)(CD,ENO) schemes + 3 ENO-LS method. 
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40 F 



- 3 2D (ENO-LS) Jacobi steps - 


Figure 14: (3,2)CD scheme + 2D ENO-LS method. 
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